High-flux neutron generation by laser-accelerated ions from single- and double-layer targets

Contemporary ultraintense, short-pulse laser systems provide extremely compact setups for the production of high-flux neutron beams, such as those required for nondestructive probing of dense matter, research on neutron-induced damage in fusion devices or laboratory astrophysics studies. Here, by coupling particle-in-cell and Monte Carlo numerical simulations, we examine possible strategies to optimise neutron sources from ion-induced nuclear reactions using 1-PW, 20-fs-class laser systems. To improve the ion acceleration, the laser-irradiated targets are chosen to be ultrathin solid foils, either standing alone or preceded by a plasma layer of near-critical density to enhance the laser focusing. We compare the performance of these single- and double-layer targets, and determine their optimum parameters in terms of energy and angular spectra of the accelerated ions. These are then sent into a converter to generate neutrons via nuclear reactions on beryllium and lead nuclei. Overall, we identify configurations that result in neutron yields as high as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^{10}\,{\mathrm{n}}\,{\mathrm{sr}}^{-1}$$\end{document}∼1010nsr-1 in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 1$$\end{document}∼1-cm-thick converters or instantaneous neutron fluxes above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{23}\,{\mathrm{n}}\,{\mathrm{cm}}^{-2}\,{\mathrm{s}}^{-1}$$\end{document}1023ncm-2s-1 at the backside of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim 100$$\end{document}≲100-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu$$\end{document}μm-thick converters. Considering a realistic repetition rate of one laser shot per minute, the corresponding time-averaged neutron yields are predicted to reach values (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gtrsim 10^7\,{\mathrm{n}} \,{\mathrm{sr}}^{-1}\,{\mathrm{s}}^{-1}$$\end{document}≳107nsr-1s-1) well above the current experimental record, and this even with a mere thin foil as a primary target. A further increase in the time-averaged yield up to above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^8\,{\mathrm{sr}}^{-1}\,{\mathrm{s}}^{-1}$$\end{document}108sr-1s-1 is foreseen using double-layer targets.


Results
Simulation setup. Our methodology is outlined in Fig. 1. We first use the particle-in-cell (PIC) calder code (see "Methods") to simulate in 2D3V (two-dimensional in space, three-dimensional in momentum space) geometry, the laser-matter interaction and proton (or deuteron) acceleration from either SLTs or DLTs. In a second step, the accelerated ions are transferred to a three-dimensional (3D) Monte Carlo code which describes neutron generation in a secondary converter target (see "Methods").
The laser parameters in the PIC simulations are chosen to match those already accessible at the Apollon facility 24 , but assuming improved temporal contrast conditions-such as those expected from fielding a plasma mirror system 29 -in order to enable efficient interaction of the laser pulse with nanometer-scale foils, as investigated in the following. The laser is modeled as a Gaussian pulse of central wavelength L = 0.8 µm and τ L = 20 fs full-width-at-half-maximum (FWHM) duration. It is focused to a D L = 5 µm FWHM spot at the front side of the target (i.e. the front side of the NCD plasma layer in the case of a DLT). Its peak intensity www.nature.com/scientificreports/ is I 0 = 2 × 10 21 W cm −2 , corresponding to a dimensionless field strength a 0 = eE 0 /m e cω L = 30.6 ( E 0 is the laser electric field, c the light speed, e the elementary charge, m e the electron mass, and ω L = 2πc/ L the laser frequency). PIC simulations are performed in the x − y plane. The laser pulse is linearly polarised along y and propagates in the +x direction. In a realistic 3D geometry, the laser pulse energy and power would be of 22 J and 1 PW, close to the current Apollon parameters. As detailed in Table 1, the simulated DLTs consist of a submicron-thick, fully ionised plastic CH 2 (or CD 2 ) foil of solid density, preceded by a fully ionised carbon (C 6+ ) NCD plasma layer of varying density ( n e,NCD ) and length ( l NCD ) . The NCD parameters optimizing proton acceleration from DLTs have been investigated both numerically and analytically by Pazzaglia et al. 43 . In that study, the maximum proton cutoff energies were attained for NCD layers of thickness close to the relativistic self-focusing length, and over a limited density range. The following approximate formulas were obtained for the optimal length and density of the NCD layer [see Eqs. (23) and (24) in Ref. 43 ]: Figure 2. Relativistic self-focusing of a 2 × 10 21 W cm -2 , 20 fs duration laser pulse in near-critical plasmas. Spatial distributions of (a) the electron density n e (saturated colormap) and (b) E y electric field component for an initial plasma electron density n e,NCD = 1.06 n cr at the simulation time t = 136.9 fs , i.e., 99.3 fs after the laser pulse maximum entered the simulation box. (c) Time evolution of the instantaneous peak intensity of the laser pulse during its propagation through the plasma for eight different values of n e,NCD . The horizontal dashed line indicates the peak laser intensity in vacuum, I 0 = 2 × 10 21 W cm −2 . (d) Dependence of the maximum intensity increase factor I max /I 0 (blue triangles and solid line), focusing distance l f (green circles and dashed line), and FWHM spot size D m (red crosses and solid line) of the laser pulse on the initial plasma density. The symbols represent the PIC simulation results while the curves correspond to the formulas given by Pazzaglia et al. 43 . The vertical dashed line indicates the plasma density yielding the best focusing according to analytical theory. Table 1. Parameters of the proton acceleration simulations. The physical parameters of the single-and doublelayer simulations are: Thickness of the near-critical density (NCD) carbon plasma layer ( l NCD ), thickness of the solid-density (SD) CH 2 (or CD 2 ) layer ( l SD ), electron density of the NCD layer ( n e,NCD ), electron density of the SD layer n e,SD = 200 n cr , laser spot size at the front of the NCD layer (if present, otherwise at the front of the SD layer) D L = 5 µm , dimensionless laser field strength a 0 = 30.6 , total number of accelerated protons or deuterons with energy above 1 MeV ( N p , N d ), ion cutoff energy ( E cutoff ).  www.nature.com/scientificreports/ where γ 0 = 1 + a 2 0 /2 is the mean Lorentz factor of the laser-driven electrons and n cr [cm −3 is the nonrelativistic electron critical density. For our parameters, γ 0 = 21.7 , n e,NCD = 1.93 n cr and l NCD = 14.1 µm . These values will provide reference conditions for our PIC simulations.
Laser self-focusing in the near-critical plasma. To start with, we inspect the process of laser self-focusing in the NCD plasma layer. To illustrate it, we consider the case of n e,NCD = 1.06n c . Figure 2a,b depict the spatial distributions of, respectively, the electron density ( n e ) and transverse electric field ( E y ) at a time ( t = 151 fs ) when the laser has propagated about 21 µm in the plasma. An electron density channel has then formed along the laser path: the laser pulse is concentrated into a D m = 1.1 µm FWHM spot at the head of the channel, where its peak intensity reaches ∼ 7.4 × 10 21 W cm −2 -an increase by a factor of ∼ 3.7 over its incident value. Inside the channel, we find (not shown) that a large fraction of the plasma electrons are accelerated by the laser wave and/or wakefields 47,48,53,58 to Lorentz factors as high as γ ≃ 300-400, largely exceeding the standard ponderomotive scaling, �γ � ≃ γ 0 . Figure 2c plots the temporal evolution of the maximum intensity I max = cε 0 2 (E 2 x + E 2 y + E 2 z ) achieved within the simulation box for eight different initial electron densities in the NCD plasma, in the range 0.45 ≤ n e,NCD /n cr ≤ 2.10 . All cases lead to a significant (by a factor of 3 ), and quite comparable, enhancement of the intensity compared to its value in vacuum ( I 0 = 2 × 10 21 W cm −2 , shown as a dashed line). The maximum intensity is reached after an interaction time of ∼ 120-180 fs, decreasing with density. The n e,NCD = 1.06 n cr case, corresponding to Fig. 2a,b, is found to yield, though by only a very small margin, the highest intensity amplification. Figure 2d shows the variations with initial plasma density n e,NCD in the laser intensity amplification factor I max /I 0 as well as in the associated focal spot D m and length l f . For each value of n e,NCD , those quantities are recorded at the time and location corresponding to the peak instantaneous laser intensity. For comparison, the theoretical estimates of D m and l f , based on the thin-lens approximation, and given by Eqs. (2) and (4) of Ref. 43 are displayed as solid red and dashed green curves, respectively. Also overlaid (blue solid line) is the intensity enhancement factor I max /I 0 , as obtained by solving numerically Eqs. (7) and (8) of Ref. 43 . The dashed line indicates the density ( n e, NCD = 1.74 n cr ) that maximises the intensity enhancement ( I max /I 0 = 4.1 ) according to that model.
Overall, the stronger lensing effect of the plasma at larger density is clearly demonstrated, and the theoretical predictions match well with the simulation data. Consistent with Fig. 2c, the simulated intensity enhancement is found to weakly vary (by ∼ 3.0-3.7) in the density range investigated, 0.45 ≤ n e,NCD /n cr ≤ 2.1 . While I max /I 0 slightly decreases (from ∼ 3.7 to ∼ 3.3 ) when the plasma density is raised from n e,NCD = 1.06 n cr to 2.1 n cr , it drops relatively more abruptly when the plasma density is decreased from the optimal density of 1.06 n cr (down to I max /I 0 ≃ 3.0 at n e,NCD = 0.45).
It is worth noting that the observed variation in I max is weaker than that one would naively infer assuming laser energy conservation (and hence I max ∝ D −2 m ) from the concomitant variation in spot size. As the latter quantity shrinks from D m ≃ 1.5 µm to ≃ 0.8 µm as n e,NCD increases, one would then expect I max to increase by a factor ∼ (1.8/0.8) 2 ∼ 3.5 , at odds with the simulation data. This discrepancy points to the strong dissipation undergone by the laser pulse during its propagation in the NCD plasma. Actually, in their model Pazzaglia et al. made an attempt at taking into account laser dissipation into hot electron generation. While the simulations and the model yield comparable maximum values of the intensity increase ( I max /I 0 ≃ 3.7 vs ≃ 4.1 ), the corresponding optimal plasma densities are significantly lower in the simulations ( n e,NCD = 1.06 n cr ) than theoretically predicted ( n e,NCD ≃ 1.5 − 2 n cr ). These results suggest that, for our simulation parameters, the actual dissipation is higher than described analytically, possibly due to fast electrons being energized much beyond the ponderomotive level considered in the model 43 . For instance, the excitation of strong plasma waves could provide an efficient additional source of energy depletion of the laser pulse 48,61 .
We conclude this part by mentioning that in an actual 3D configuration, relativistic laser self-focusing is expected to lead to even stronger intensification of the laser pulse, as previously shown in Ref. 43 .
Proton acceleration from single-and double-layer targets. We now assess the performance of SLTor DLT-based ion acceleration setups in generating intense neutron fluxes from secondary converter targets. Table 1 summarises the parameters of the six SLT and DLT simulations we have carried out. Runs #1 and #2 represent SLT configurations. The thickness of the CH 2 foil in run #1 is that predicted to optimise the proton cutoff energy under our irradiation conditions according to Refs. 34,62 [see Eq. (4) in "Methods"]. The 115-nmthick CH 2 foil in run #2 is optimised for the maximum intensity of 7.8 × 10 21 W cm −2 achieved during the laser propagation in the NCD plasma. This foil is about twice thicker than in run #1 but still partly transparent; a higher number of protons should then be produced but at lower energies. Runs #3-5 use DLTs with NCD layers of different densities or lengths, coated on the same CH 2 foil as in run #2 to account for plasma lensing. In runs #3 and #4, the NCD layer thickness is set to the laser focusing length in order to achieve the highest possible www.nature.com/scientificreports/ laser intensity, and thus to maximise proton energies. The combinations of the n e,NCD and L NCD values in runs #3 and #4 are those yielding the strongest laser intensities [see Fig. 2c]. Run #5 is similar to run #3 but uses a shorter L NCD to examine how the final proton and neutron beams depend on the NCD layer length. Finally, run #3b has the same parameters as run #3 but protons are here fully replaced with deuterons (i.e. the solid foil is made of CD 2 ). We acknowledge that our comparison of DLTs with CH 2 and CD 2 solid foils is somewhat idealised because, in reality, (i) there would be proton-rich contaminant layers on either side of the foil, whatever its composition and (ii) deuterated targets used for laser-plasma acceleration often contain a 1% fraction of hydrogen within their bulk 16 . In all PIC simulations, a virtual detector is placed at a distance of 26.2 µm behind the rear side of the target. When an accelerated proton crosses this "plane", its position, momentum and time of arrival are stored. A relatively close detector position is chosen in order to interrupt the acceleration process. This is a common practice because 2D PIC simulations are known to overestimate the final proton energy 63 . Our choice of distance relies on a comparison between 2D and 3D PIC simulations 34 . Figure 3 displays the spatial distributions of the proton number density (a,c) and average energy (b,d), 207 fs after the laser pulse maximum has hit the solid foil. The top and bottom rows correspond to SLT run #2 and DLT run #4, respectively. In both cases, the target protons are fully evacuated from the laser spot region and preferentially accelerated in the forward direction due to the combined effects of RPA and TNSA (see Fig. 3a,c). Compared to the SLT case, the tighter laser focusing achieved in the DLT translates into a more energetic (by a factor of ∼ 3 ), yet more divergent, proton beam originating from a narrower source (cf. Fig. 3b,d). Figure 4a compares the energy spectra of the protons (or deuterons) that reached the detector plane [marked by vertical dashed lines in Fig. 3 for the six configurations considered. The absolute maximum proton energy ( ∼ 240 MeV ) is recorded in DLT runs #3 and #4, as expected from the associated laser intensification (see Fig. 2). By comparison, the SLT cases yield significantly lower proton energies: ∼ 130 MeV in the optimised run #1 and ∼ 95 MeV for the about twice thicker foil of run #2.  www.nature.com/scientificreports/ As already mentioned, the flip side of the fastest protons produced in DLTs is their increased angular spread, due to strong focusing and partial transmission of the laser pulse through the foil. Actually, their angular spectra plotted in Fig. 4b, turn out to be double-peaked at ∼ ± 15-20 • and to extend up to angles as large as ∼ ± 40 • . By contrast, SLTs can provide more collimated protons and also in greater numbers (see Table). Run #2 indeed yields a proton angular spectrum peaking on axis and with an FWHM spread of ∼ 12.5 • .
Deuterons in DLT run #3b attain a similar, albeit slightly lower, cutoff energy ( ∼ 225 MeV ) than protons in a similar setup (run #3), but their number is about 2-4× lower in the high-energy ( > 100 MeV ) part of the spectrum. Their angular distribution also peaks off axis but with a less pronounced on-axis minimum than in the proton cases. We note that while the absence of surface protons in our simulations is expected to favour, to some extent, the deuteron acceleration, it was shown experimentally 16 that a large majority of the ions accelerated from 300 nm CD 2 foils were deuterons despite a deuterisation level only as high as 90%.
Proton transport and neutron generation through the converter. The accelerated protons (or deuterons) recorded by the virtual detector in the PIC simulations are used as inputs to the Monte Carlo simulations of their transport through the neutron converter target. The latter, located 26.2 µm away from the laser target, consists of a lead ( 208 Pb) or beryllium ( 9 Be) cylinder of 3-cm radius and varying length (l). This procedure is illustrated in Fig. 5. Figure 5a displays the time-resolved energy spectrum of the outgoing protons in SLT case #2: it is characterised by a maximum energy of ≃ 95 MeV (i.e., the lowest cutoff energy among our simulations), and a root-mean-square pulse duration of ≃ 190 fs (the time interval between the incidence of the first and last protons with energy above 1 MeV is 770 fs). The dashed orange curve in Fig. 5b plots the fraction of energy lost (via inelastic collisions) by those protons in the Pb converter, as a function of its length. The dashed blue curve plots the same quantity for the proton beam generated in DLT case #4, associated with the highest ( ∼ 250 MeV ) cutoff energy. The orange and blue solid lines represent, respectively for cases #2 and #4, the (complementary) fraction of proton energy transmitted across the converter's backside. For case #2 (resp. case #4), the beam energy dissipation remains negligible ( ≤ 1%) in Pb converters thinner than ≃ 250 µm (resp. ≃ 2 mm ), while it is almost complete ( ≥ 90% ) for l ≥ 3 mm (resp. l ≥ 3 cm). Figure 6 details the variations with the converter length of the properties of the neutron sources generated in beryllium (left column) and lead (right column), for the ion-acceleration setups listed in Table 1. Only the neutrons leaving the backside of the converter are recorded (i.e., those escaping from the front and lateral sides are excluded). Figure 6a,b show that the highest absolute neutron yield, N n ≃ 2.8 × 10 10 , is achieved in a 1-cmthick Pb converter exposed to the protons generated in DLT run #3-the setup generating the greatest number of fast protons. Other DLT runs #4 and #5 give very close results. The higher energy and number of protons accelerated in DLTs translate into a ∼ 40% greater total yield in Pb than in Be. The enhanced performance of Pb converters coupled with DLTs can be partly ascribed to a larger 208 Pb(p, n) cross section, which reaches ∼ 1 b at a ∼ 13 MeV proton energy and keeps on rising at higher energy, exceeding ∼ 10 b at ∼ 100 MeV 64 . This is unlike the 9 Be(p, n) cross section which peaks around ∼ 5 MeV with a value of ∼ 0.2 b , and drops at higher energy, falling below ∼ 10 mb above ∼ 100 MeV 65 .
When using SLTs, by contrast, the total neutron yield is in general larger in Be, by tens of per cent in run #1 and by up to three times in run #2 the one giving the slowest protons). The maximum yield ( N n ≃ 1.2 × 10 10 ) is then recorded in a 0.8-cm-thick Be converter using proton source #1.
For a given ion beam, the N n vs. l curve peaks around the beam penetration distance in the converter. The decreasing trend at larger l mainly originates from neutrons escaping from the lateral sides of the converter; such side losses become significant when the transverse size of the neutron beam, w n (plotted vs. l in Fig. 6c,d), approaches the converter's diameter. The beam size is obtained by fitting to a Gaussian the neutron fluence profile at the backside of the converter. Within the range of converter lengths considered, w n is approximately equal to l, except in very thin converters ( l 100 µm ), where it is mainly set by the initial transverse proton beam size and divergence to a lower value of ∼ 150 µm . The impact of the converter's shape and size on the neutron yield in accelerator-based spallation neutron sources is discussed in Ref. 31 .  Table 1). (b) Dashed lines: dissipated fraction of proton beam energy due to inelastic collisions as a function of the Pb converter length l, in SLT case #2 (orange) and DLT case #4 (red). Solid lines: transmitted fraction of proton beam energy in SLT case #2 (orange) and DLT case #4 (red). www.nature.com/scientificreports/ In Pb converters of thickness l < 1 cm , the total neutron yield is maximised in DLT runs #3 and #5, while the three DLT setups give similar results for l > 1 cm . For thin Be converters ( l 2 mm ), the neutron yield is maximised, to within small differences, in SLT runs #1 and #2 and DLT runs #3 and #5. For 0.2 < l 3 cm , runs #3 and #5 perform the best and produce the maximum absolute value N n ≃ 2 × 10 10 around l ≃ 2 cm . In thicker Be converters N n steadily decreases and is optimised in the deuteron-based run #3b.
The peak neutron flux, F n , plotted in Figs. 6e,f is a relevant parameter for certain purposes such as laboratory studies on r−process nucleosynthesis 6,66,67 . It is calculated as where dN n /dt is obtained from our extended fluka routine at the converter's backside. It is seen that neutron fluxes exceeding 10 23 n cm −2 s −1 can be attained using Be or Pb converters a few 10 µm thick only, in which the neutron source is the narrowest and the shortest. In an actual experiment, it would then be important to place the converter as close as possible to the ion-generating target, while not hindering the ion acceleration process, as is the case in our simulation setup. Note that target assemblies consisting of two solid foils separated by tens The record peak flux ( ∼ 6 × 10 23 n cm −2 s −1 ) is achieved, though by a very short margin, in DLT case #3 with a ∼ 25-µm-thick Pb converter. Other DLT cases yield very similar results while the SLT cases perform only slightly less well. For ∼ 1 cm converter lengths maximizing the neutron yield, the peak flux drops to 10 18 − 10 19 n cm −1 s −1 , with a larger difference between SLT and DLT data in Pb than in Be.
To assess the dependence of the neutron flux on the distance (d) between the laser target and the converter, we have performed two additional Monte Carlo simulations using the ion source from DLT case #3 but with a Pb converter located 1 mm away from the laser target. The peak flux at the backside of a thin ( 25 µm ) converter is then reduced to F n ≃ 9.1 × 10 22 n cm −2 s −1 , i.e., ∼ 15 % of the value achieved in our baseline configuration ( d = 26 µm ). For a 2-cm-thick converter, we obtain F n ≃ 2.1 × 10 18 n cm −2 s −1 , i.e., ∼ 75 % of the value attained previously. Figure 6g,h display the neutron yield per unit solid angle, n , i.e., the important quantity for most applications. n is measured on axis 20 cm away from the rear side of the converter. When using protons, n is maximised with DLTs. Specifically, one finds that n reaches similar peak values ( ≃ 8 × 10 9 n sr −1 ) in ∼ 2-3-cm-thick Be or Pb converters exposed to proton beams from DLT case #3. It is worth noting that the lower-energy protons from both SLT cases are only slightly less efficient in Be ( n ≃ 2.5-4 ×10 9 n sr −1 ) than DLT protons but perform significantly poorer in Pb ( n ≃ 0.8-2 ×10 9 n sr −1 ). A major finding, however, is that for Apollon-class laser parameters, the highest absolute neutron yield per solid angle ( n ≃ 1.6 × 10 10 n sr −1 ) is achieved by deuterons from DLT run #3b in a 3-cm-thick Be converter (see Fig. 6g). This twofold increase in n may seem surprising as deuterons produce an approximately three times lower neutron yield than DLT-accelerated protons in Pb. The origin of this result is the pronounced forward directionality of deuteron breakup neutrons for deuteron energies above ∼ 2.2 MeV 69,70 . Figure 7a depicts the neutron fluence profile through the converter when coupling proton source #3 with a 2-cm-thick Pb converter (i.e. the proton-based setup maximizing N n and n in Pb). The front and rear boundaries of the converter are indicated by white dashed lines. Neutron fluences exceeding 10 12 n cm −2 are found at depths x 2 mm around the proton beam axis. At the backside ( x = 2 cm ), the neutron fluence drops to ∼ 10 10 n cm −2 over a ∼ 2.3 cm FWHM transverse width. Figure 7b shows the corresponding energy-angle neutron spectrum. As only neutrons escaping from the rear side of the converter are considered, the angular spectrum is restricted to solid angles ≤ 2π . Two neutron populations can be distinguished: (i) an approximately isotropic group of low-energy ( 3 MeV ) neutrons, which amounts to ∼ 70% of all generated neutrons; (ii) a group of much more energetic neutrons (up to ∼ 160 MeV ), the divergence of which decreases at higher energy. Figure 8 details the temporal evolution of the neutron distribution produced from a 100-µm-thick Pb converter by proton source #3. This configuration generates a total number of N n ≃ 1.2 × 10 9 neutrons. The plots in the left-and right-hand side columns visualise, respectively, the short ( ≤ 50 ps ) and long ( ≤ 3 ns ) duration profiles of the neutron source. The early-time neutron burst, with a peak instantaneous flux of F n ≃ 2.3 × 10 23 n cm −2 s −1 , contains the most energetic neutrons (up to ∼ 230 MeV ) and is delivered within a few picoseconds (see Fig. 8a,b). The outgoing neutron flux steadily decreases afterwards, yet remains above ∼ 10 20 n cm −2 s −1 (resp. ∼ 10 19 n cm −2 s −1 ) till t ≃ 270 ps (resp. t ≃ 0.9 ns).

Discussion
By coupling particle-in-cell and Monte Carlo simulations, we have investigated numerically the possibility of generating high-flux neutron sources by 1 PW-class, 20 fs laser pulses as are now available at the Apollon facility 24 . Such sources rely on nuclear reactions triggered in a light ( 9 Be) or heavy ( 208 Pb) converter by fast (up to ∼ 100-200MeV ) ions (protons or deuterons) driven from the laser target. Inspired by current trends in the laser-plasma community, we have examined the potential of double-layer targets-comprising a few-micronscale, near-critical plasma layer attached to a nanometric-scale solid foil-in enhancing ion acceleration and  Table 1). The yellow line is a Gaussian fit (with a FWHM width w n = 2.3 cm ) of the neutron fluence as recorded at the converter's backside ( x = 2 cm ). (b) Energy-angle spectrum of the outgoing neutrons. www.nature.com/scientificreports/ the resulting neutron production. Our simulations predict that in specially designed DLTs, the peak intensity of the plasma-focused laser pulse can be increased nearly fourfold, entailing more than doubled maximum proton energies compared to those obtained with plain thin solid foils. This translates into approximately fourfold increased neutron yields (whether measured per unit solid angle or angle-integrated) from cm-scale Pb converters, in which ∼ 10 − 100 MeV protons benefit from large ( ∼ 1-10b ) (p, n) cross sections. It should be noticed that besides intensifying the laser light and boosting proton acceleration, the NCD layer in DLTs can be beneficial in shielding the ultrathin solid foil against the laser pedestal or prepulses 71 . By contrast, when employing a Be converter, more comparable neutron yields are predicted to be released by proton beams accelerated from SLTs and DLTs, due to decreasing (p, n) cross sections at proton energies > 5 MeV . Interestingly, we find that the overall maximum neutron yield per unit solid angle achieved with protons (be it in Pb or Be converters) can be surpassed about twofold with deuterons by using a deuterated DLT and a Be converter.
Our simulation study also indicates that Apollon-class systems should be capable of generating peak neutron fluxes in excess of 10 23 n cm −2 s −1 using either SLTs or DLTs. Converter targets of a few tens of microns would then be required to limit the pulse duration and lateral spread of the emitted neutrons.
The above sources, and especially those utilising DLTs, compare favourably with previously reported experimental [16][17][18] or numerical 23 works on laser-driven neutron sources but using much longer ( τ L ≃ 0.5-1ps ) and more energetic ( ∼ 4-14× ) laser pulses. This is shown in Table 2 which summarises the results of these past works and ours, detailing in each case the laser parameters, the nature of the projectile and converter material, and the (measured or simulated) features of the neutron source.
The highest neutron yield per unit solid angle obtained in our study ( n ≃ 1.6 × 10 10 n sr −1 in deuteronbased DLT case #3b) is very close to the current experimental record ( n ≃ 1.4 × 10 10 n sr −1 ) reported in Ref. 17 . Moreover, our optimum proton-and deuteron-based setups #3 and #3b are predicted to produce neutron yields per unit laser pulse energy of � n /E L ≃ 3.4 − 7.0 × 10 8 n sr −1 J −1 , exceeding the record-high values ( � n /E L ≃ 0.95-2.5 ×10 8 n sr −1 J −1 ) achieved at the PHELIX facility 17,18 , and performing similarly to the scheme (based on radiation pressure acceleration in overcritical CD 2 foams by a circularly polarised laser pulse) recently proposed in Ref. 23 .
Our numerical results appear even more promising as regards the time-averaged neutron yield per unit solid angle, that is, the product of the neutron yield and the laser shot repetition rate. This comparison evidently implies that our target setups can be adapted to the relatively high repetition rate (one shot per minute) allowed by Apollon-class systems. When combining the CH 2 SLT #1 and a Pb converter of optimum length ( ∼ 1 cm ), a time-averaged neutron yield of f L n ≃ 3 × 10 7 n sr −1 s −1 is expected, a value about ten times higher than measured experimentally so far. When attaching an optimum-density plasma layer to the optimum-thickness solid CH 2 layer (DLT case #3), this value can rise to f L n ≃ 1.2 × 10 8 n sr −1 s −1 . Finally, using the same DLT Simulations of proton acceleration. Ion acceleration from both SLTs and DLTs was modeled in 2D3V geometry with the PIC calder code. The size of the simulation domain was L x × L y = 63.7 × 50.9 µm 2 , discretised into 10,000 ×8000 cells with x = y = 6.37 nm . A full 3D simulation of the problem with the same level of discretisation remains well outside our computational reach. The laser pulse maximum entered the simulation box at t = 40.32 fs . These simulations captured both the self-focusing of the laser pulse through the NCD carbon plasma layer (if present) and its interaction with the solid-density (SD) CH 2 or CD 2 layer, from which originate the protons or deuterons used for neutron production.
The SD layer was modeled by 50 macroparticles per cell and species (electrons, C 6+ , and H + or D + ) while 15 macroparticles per cell and species were used in the NCD layer (electrons, C 6+ ). Fourth-order shape factors were used for the macroparticles.
According to the PIC simulation results of Refs. 34,62 , the thickness of the SD layer that optimises the cutoff ion energy with a femtosecond laser pulse is given by where n e,SD is the electron density of the SD layer. For the density n e,SD = 200 n cr used in our simulations, we obtain l opt ≃ 62 nm . This value served as a reference to design our SD targets.
Simulations of neutron generation. The 3D Monte Carlo fluka code 73,74 was employed to describe the neutron generation from the fast protons reaching the virtual detector in the 2D3V calder PIC simulation. To this goal, we extended fluka with two new modules allowing the code to accept as inputs the macro-particles tracked in the PIC simulation.
The first module enables the input particles to be injected with the same temporal profile as recorded in the PIC simulation. The second module converts the set of macro-particles into an axisymmetric distribution. In detail, each macro-particle's position (x, y) and momentum (p x , p y ) is rotated by a random azimuthal angle. Moreover, to obtain the number of physical particles represented by the macro-particle in 3D geometry, its original (2D) statistical weight (a linear density in 2D geometry) is multiplied by 2πy 0 , where y 0 is its initial distance to the axis in the unperturbed target (prior to the laser irradiation). A similar procedure was adopted by Jiang et al. 75 to inject a 2D-PIC-simulated electron distribution into a 3D hybrid-PIC code. We note that this post-processing modifies the absolute laser-to-ion conversion efficiency: it is evaluated to be of ∼ 9% which is higher than, yet comparable with, the value of ∼ 5% directly extracted from the 2D simulation. These values are (4) l opt ≃ 0.5a 0 n cr n e,SD L , Table 2. Comparison of the optimum neutron sources studied in this work with state-of-the-art published results. The first column details the reference or simulation case considered (experimental works are highlighted in bold). The next three columns present the corresponding laser pulse energy ( E L ), duration ( τ L ), and the shot frequency ( f L ). The fifth column details the ion projectile and the converter material, and the last three columns summarise the total neutron yield per unit solid angle ( n ), the neutron yield normalised to the laser energy ( � n /E L ) and the time-averaged neutron yield ( f L n ).
Reference or simulation # E L (J) τ L (fs) f L (Hz) Configuration n ( n sr −1 ) � n /E L ( n J −1 sr −1 ) f L n ( n s −1 sr −1 ) Roth et al. 16 76 in relativistically transparent foils, albeit using a higher-energy ( ∼ 200 J ), longer-duration ( ∼ 0.9 ps ) laser drive. Each fluka simulation made use of a minimum of 3.6 × 10 8 particles to model the incident proton or deuteron flux. According to convergence tests, this choice ensured the statistical accuracy of the output data. Additional simulations were also conducted to assess the contribution to neutron production of photonuclear reactions induced by fast electrons escaping the laser target. Under our conditions, this mechanism was found to be negligible compared to proton-induced neutron generation.
As fluka does not yet include models for deuteron transport, we resorted to the Monte Carlo mcnp-6 code 77 to simulate deuteron-induced neutron generation. These simulations did not provide output on the temporal profile of the neutron pulse.
In all Monte Carlo simulations, the converter target was a 3-cm-radius, finite-length cylinder, made of beryllium or lead. The former material is that used in state-of-the-art laser experiments while the latter is representative of high-Z converters, susceptible to spallation-type neutron generation for proton energies above ∼ 200 MeV 31,78 .

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper. Simulation results are available from the corresponding author upon reasonable request.